[Group_13]_Project_CK_03

Project 3 - Body performance Data

Hiện nay các phong trào tập thể thao đang ngày một phát triển, thu hút nhiều nhóm tuổi và giới tính. Dữ liệu bodyPerformance.csv chứa thông tin của 13,393 người tham gia tập thể thao tại Hàn Quốc, với 12 biến như sau:

• age - độ tuổi (từ 20 tới 64);

• gender - giới tính (F: nữ, M: nam);

• height_cm - chiều cao (đơn vị: cm);

• weight_kg - cân nặng (đơn vị: kg);

• body fat_% - phần trăm mỡ cơ thể (%);

• diastolic - huyết áp tâm trương (phút);

• systolic - huyết áp tâm thu (phút);

• gripForce - lực kẹp;

• sit and bend forward_cm - ngồi và gập người về phía trước;

• sit-ups counts - số lần gập bụng;

• broad jump_cm - nhảy xa (đơn vị: cm);

• class - phân lớp hiệu suất (A: tốt nhất, B,C,D).

Hãy xử lý dữ liệu này để giúp cho các chuyên gia sức khỏe biết được hiệu quả của việc tập thể dục, và các yếu tố ảnh hưởng tới hiệu quả.

• Load dữ liệu

# Khai báo các thư viện
library(tidyverse)

## ── Attaching core tidyverse packages ──────────────────────── tidyverse 2.0.0 ──
## ✔ dplyr     1.1.4     ✔ readr     2.1.5
## ✔ forcats   1.0.0     ✔ stringr   1.5.1
## ✔ ggplot2   3.5.1     ✔ tibble    3.2.1
## ✔ lubridate 1.9.3     ✔ tidyr     1.3.1
## ✔ purrr     1.0.2     
## ── Conflicts ────────────────────────────────────────── tidyverse_conflicts() ──
## ✖ dplyr::filter() masks stats::filter()
## ✖ dplyr::lag()    masks stats::lag()
## ℹ Use the conflicted package (<http://conflicted.r-lib.org/>) to force all conflicts to become errors

library(readr)

# Đọc file csv
data <- read.csv(file = "D:/bodyPerformance.csv")
data <- data |> janitor::clean_names() # Chuyển tên các cột về chữ thường
glimpse(data)

## Rows: 13,393
## Columns: 12
## $ age                     <dbl> 27, 25, 31, 32, 28, 36, 42, 33, 54, 28, 42, 57…
## $ gender                  <chr> "M", "M", "M", "M", "M", "F", "F", "M", "M", "…
## $ height_cm               <dbl> 172.3, 165.0, 179.6, 174.5, 173.8, 165.4, 164.…
## $ weight_kg               <dbl> 75.24, 55.80, 78.00, 71.10, 67.70, 55.40, 63.7…
## $ body_fat                <dbl> 21.3, 15.7, 20.1, 18.4, 17.1, 22.0, 32.2, 36.9…
## $ diastolic               <dbl> 80, 77, 92, 76, 70, 64, 72, 84, 85, 81, 63, 69…
## $ systolic                <dbl> 130, 126, 152, 147, 127, 119, 135, 137, 165, 1…
## $ grip_force              <dbl> 54.9, 36.4, 44.8, 41.4, 43.5, 23.8, 22.7, 45.9…
## $ sit_and_bend_forward_cm <dbl> 18.4, 16.3, 12.0, 15.2, 27.1, 21.0, 0.8, 12.3,…
## $ sit_ups_counts          <dbl> 60, 53, 49, 53, 45, 27, 18, 42, 34, 55, 68, 0,…
## $ broad_jump_cm           <dbl> 217, 229, 181, 219, 217, 153, 146, 234, 148, 2…
## $ class                   <chr> "C", "A", "C", "B", "B", "B", "D", "B", "C", "…

• Bảng tóm tắt và khái quát về dữ liệu

library(dplyr)
# Tính toán thống kê tóm tắt
summary_table <- data %>%
  reframe(
    Bien = c("age", "height_cm", "weight_kg", "body_fat", "diastolic", "systolic", "grip_force", "sit_and_bend_forward_cm", "sit_ups_counts", "broad_jump_cm"),
    n = n(),
    Trung_binh = c(mean(age), mean(height_cm), mean(weight_kg), mean(body_fat), mean(diastolic), mean(systolic), 
                   mean(grip_force), mean(sit_and_bend_forward_cm), mean(sit_ups_counts),   mean(broad_jump_cm)),
    Trung_vi = c(median(age), median(height_cm), median(weight_kg), median(body_fat), median(diastolic), median(systolic),
                 median(grip_force), median(sit_and_bend_forward_cm), median(sit_ups_counts), median(broad_jump_cm)),
    Min = c(min(age), min(height_cm), min(weight_kg), min(body_fat), min(diastolic), min(systolic),
            min(grip_force), min(sit_and_bend_forward_cm), min(sit_ups_counts), min(broad_jump_cm)),
    Max = c(max(age), max(height_cm), max(weight_kg), max(body_fat), max(diastolic), max(systolic),
            max(grip_force), max(sit_and_bend_forward_cm), max(sit_ups_counts), max(broad_jump_cm)),
  )

# Hiển thị bảng
library(knitr)
kable(summary_table, col.names = c("Biến", "Số lượng", "Trung bình", "Trung vị", "Min", "Max"), align = "c")

Biến	Số lượng	Trung bình	Trung vị	Min	Max
age	13393	36.77511	32.0	21.0	64.0
height_cm	13393	168.55981	169.2	125.0	193.8
weight_kg	13393	67.44732	67.4	26.3	138.1
body_fat	13393	23.24017	22.8	3.0	78.4
diastolic	13393	78.79684	79.0	0.0	156.2
systolic	13393	130.23482	130.0	0.0	201.0
grip_force	13393	36.96388	37.9	0.0	70.5
sit_and_bend_forward_cm	13393	15.20927	16.2	-25.0	213.0
sit_ups_counts	13393	39.77122	41.0	0.0	80.0
broad_jump_cm	13393	190.12963	193.0	0.0	303.0

• Kiểm tra và xóa các giá trị outliers

remove_outliers <- function(x) {
  Q1 <- quantile(x, 0.25)
  Q3 <- quantile(x, 0.75)
  IQR <- Q3 - Q1
  lower_bound <- Q1 - 1.5 * IQR
  upper_bound <- Q3 + 1.5 * IQR
  x[x >= lower_bound & x <= upper_bound]
}

# Remove outliers within each class
data <- data %>%
  group_by(class) %>%
  filter(weight_kg >= quantile(weight_kg, 0.25) - 1.5 * IQR(weight_kg) & 
           weight_kg <= quantile(weight_kg, 0.75) + 1.5 * IQR(weight_kg) &
           height_cm >= quantile(height_cm, 0.25) - 1.5 * IQR(height_cm) & 
           height_cm <= quantile(height_cm, 0.75) + 1.5 * IQR(height_cm) &
           body_fat >= quantile(body_fat, 0.25) - 1.5 * IQR(body_fat) &
           body_fat <= quantile(body_fat, 0.75) + 1.5 * IQR(body_fat) &
           diastolic >= quantile(diastolic, 0.25) - 1.5 * IQR(diastolic) &
           diastolic <= quantile(diastolic, 0.75) + 1.5 * IQR(diastolic) &
           systolic >= quantile(systolic, 0.25) - 1.5 * IQR(systolic) &
           systolic <= quantile(systolic, 0.75) +1.5*IQR(systolic) &
           grip_force >= quantile(grip_force, 0.25) -1.5*IQR(grip_force) &
           grip_force <= quantile(grip_force, 0.75) +1.5*IQR(grip_force) &
           sit_and_bend_forward_cm >= quantile(sit_and_bend_forward_cm, 0.25) - 1.5*IQR(sit_and_bend_forward_cm) &
           sit_and_bend_forward_cm <= quantile(sit_and_bend_forward_cm, 0.75) +1.5*IQR(sit_and_bend_forward_cm) &
           sit_ups_counts >= quantile(sit_ups_counts, 0.25) -1.5*IQR(sit_ups_counts) &
           sit_ups_counts <= quantile(sit_ups_counts, 0.75) +1.5*IQR(sit_ups_counts) &
           broad_jump_cm >= quantile(broad_jump_cm, 0.25) -1.5*IQR(broad_jump_cm) &
           broad_jump_cm <= quantile(broad_jump_cm, 0.75) +1.5*IQR(broad_jump_cm) 
  )

Việc lọc outliers cho toàn bộ dữ liệu mà không theo nhóm class làm cho dữ liệu của class D bị sai sót. Do đó, ta lọc outliers dựa trên dữ liệu của mỗi nhóm, việc này giúp giữ nguyên được tính chất của dữ liệu ban đầu.

• Vẽ boxplot so sánh các biến theo phân lớp hiệu suất (class)

library(ggplot2)

# Danh sách các biến cần so sánh
variables <- c("age", "height_cm", "weight_kg", "body_fat", "diastolic", "systolic", 
               "grip_force", "sit_and_bend_forward_cm", "sit_ups_counts", "broad_jump_cm")

# Tạo boxplot cho từng biến
plots <- lapply(variables, function(var) {
  ggplot(data, aes_string(x = "class", y = var, fill = "class")) +
    geom_boxplot() +
    labs(title = paste("Distribution of", var, "by Class"),
         x = "Class", y = var) +
    theme_minimal()
})

## Warning: `aes_string()` was deprecated in ggplot2 3.0.0.
## ℹ Please use tidy evaluation idioms with `aes()`.
## ℹ See also `vignette("ggplot2-in-packages")` for more information.
## This warning is displayed once every 8 hours.
## Call `lifecycle::last_lifecycle_warnings()` to see where this warning was
## generated.

# Hiển thị các biểu đồ
for (plot in plots) {
  print(plot)
}

• Bảng tóm tắt và khái quát dữ liệu phân theo lớp hiệu suất

custom_summary <- function(x) {
  data.frame(
    n = length(x),
    mean = mean(x, na.rm = TRUE),
    sd = sd(x, na.rm = TRUE),
    median = median(x, na.rm = TRUE),
    trimmed = mean(x, trim = 0.1, na.rm = TRUE),
    mad = mad(x, na.rm = TRUE),
    min = min(x, na.rm = TRUE),
    max = max(x, na.rm = TRUE),
    range = max(x, na.rm = TRUE) - min(x, na.rm = TRUE)
  )
}

grouped <- split(data[, -c(2, 12)], data$class)
result <- lapply(grouped, function(group) {
  sapply(group, custom_summary)
})

print(result)

## $A
##         age      height_cm weight_kg body_fat diastolic systolic grip_force
## n       3303     3303      3303      3303     3303      3303     3303      
## mean    35.2198  167.8805  64.40722  20.51164 77.92928  129.2894 38.62281  
## sd      12.96312 7.806794  10.50914  6.336261 10.39194  14.25649 10.87522  
## median  30       168       64.5      20.2     78        128      39.1      
## trimmed 33.68558 167.8977  64.17427  20.34585 77.96913  129.1154 38.36367  
## mad     10.3782  8.74734   12.30558  6.81996  10.3782   14.826   14.23296  
## min     21       145       34.5      3        49        88       2.1       
## max     64       190.9     97        38.8     106       167      70.5      
## range   43       45.9      62.5      35.8     57        79       68.4      
##         sit_and_bend_forward_cm sit_ups_counts broad_jump_cm
## n       3303                    3303           3303         
## mean    21.28712                47.90675       202.8707     
## sd      4.179923                10.80746       35.81801     
## median  21.1                    49             202          
## trimmed 21.1539                 48.27696       203.4362     
## mad     4.29954                 11.8608        44.478       
## min     11.8                    17             94           
## max     32.9                    80             299          
## range   21.1                    63             205          
## 
## $B
##         age      height_cm weight_kg body_fat diastolic systolic grip_force
## n       3285     3285      3285      3285     3285      3285     3285      
## mean    37.02983 168.6511  66.62328  21.96273 78.73029  130.6545 37.97882  
## sd      13.65834 8.012087  10.68315  6.599006 10.32314  14.17329 10.3042   
## median  32       169.3     67.2      21.6     79        130      39.6      
## trimmed 35.90529 168.8185  66.52599  21.78668 78.78623  130.658  37.94992  
## mad     13.3434  8.30256   11.56428  6.96822  10.3782   14.826   12.30558  
## min     21       146.3     38.1      4.7      49        90       4.4       
## max     64       191.3     98.3      40.6     103       167      69.9      
## range   43       45        60.2      35.9     54        77       65.5      
##         sit_and_bend_forward_cm sit_ups_counts broad_jump_cm
## n       3285                    3285           3285         
## mean    17.306                  42.71702       195.6451     
## sd      4.697493                11.70452       36.10045     
## median  17                      44             200          
## trimmed 17.10103                43.09471       196.8847     
## mad     4.89258                 13.3434        41.5128      
## min     7.1                     12             84           
## max     30.5                    74             290          
## range   23.4                    62             206          
## 
## $C
##         age      height_cm weight_kg body_fat diastolic systolic grip_force
## n       3283     3283      3283      3283     3283      3283     3283      
## mean    36.63174 169.2662  66.83017  22.56194 78.57843  129.9427 36.69148  
## sd      13.73489 8.412133  10.73871  6.133272 10.35322  14.29215 10.13843  
## median  32       170.1     67.4      22.2     79        130      38.2      
## trimmed 35.42558 169.5314  66.80201  22.40848 78.69128  129.9494 36.66781  
## mad     13.3434  8.74734   11.71254  5.78214  10.3782   14.826   11.71254  
## min     21       145.5     38.1      6.5      50        91       11.3      
## max     64       191.4     97.4      38.3     107       165      65.2      
## range   43       45.9      59.3      31.8     57        74       53.9      
##         sit_and_bend_forward_cm sit_ups_counts broad_jump_cm
## n       3283                    3283           3283         
## mean    14.28909                38.85014       189.2811     
## sd      5.762414                12.65448       38.31712     
## median  13.9                    40             194          
## trimmed 14.04041                39.21888       190.8854     
## mad     6.52344                 13.3434        42.9954      
## min     2.3                     7              76           
## max     31.3                    71             303          
## range   29                      64             227          
## 
## $D
##         age      height_cm weight_kg body_fat diastolic systolic grip_force
## n       3254     3254      3254      3254     3254      3254     3254      
## mean    38.02274 168.6446  71.79461  27.61545 80.17342  131.1042 34.7941   
## sd      13.83676 8.978501  13.39088  7.284413 10.63221  14.74647 10.48771  
## median  36       169.8     71.61     27.3     80        131      35.7      
## trimmed 37.14132 168.8645  71.54618  27.58136 80.37876  131.2608 34.74869  
## mad     17.7912  9.78516   14.39605  7.413    11.8608   16.3086  12.45384  
## min     21       143.6     37.3      8.1      52        86       0         
## max     64       192       110.1     47.7     110       169      68.4      
## range   43       48.4      72.8      39.6     58        83       68.4      
##         sit_and_bend_forward_cm sit_ups_counts broad_jump_cm
## n       3254                    3254           3254         
## mean    7.770553                30.06331       174.638      
## sd      9.153477                14.92396       40.65272     
## median  7.8                     31             180          
## trimmed 7.998587                30.38364       176.3253     
## mad     9.48864                 16.3086        44.478       
## min     -17.4                   0              50           
## max     32.1                    67             275          
## range   49.5                    67             225

Nhóm A (hiệu suất cao nhất) có các đặc điểm nổi bật: Tuổi trẻ hơn, cân nặng và mỡ cơ thể thấp hơn, khả năng sức mạnh (lực kẹp, nhảy xa), sức bền (gập bụng), và độ dẻo dai cao hơn đáng kể so với các nhóm khác. Ngược lại, nhóm D (hiệu suất thấp nhất) có tuổi cao hơn, cân nặng và mỡ cơ thể lớn hơn, cùng với các chỉ số thể chất thấp hơn.

• Vẽ ma trận tương quan

data_hm <- data

# Chuyển đổi cột gender thành 0 và 1
 data_hm$gender <- ifelse(data$gender == "M", 1, 0)

# Chuyển đổi cột class thành dạng số
 data_hm$class <- as.numeric(factor(data$class, levels = c("D", "C", "B", "A")))

library(reshape2)

## 
## Attaching package: 'reshape2'

## The following object is masked from 'package:tidyr':
## 
##     smiths

cor_matrix <- cor(data_hm %>% select_if(is.numeric)) # Chọn các biến số
ggplot(melt(cor_matrix), aes(Var1, Var2, fill = value)) +
  geom_tile(color = "white") +
  scale_fill_gradient2(low = "blue", high = "red", mid = "white", midpoint = 0, 
                       limit = c(-1, 1), space = "Lab", 
                       name = "Correlation") +
  theme_minimal() +
  labs(title = "Ma trận tương quan", x = "", y = "") +
  theme(axis.text.x = element_text(angle = 45, vjust = 1, hjust = 1))

Ma trận tương quan cho thấy class có mối quan hệ tương quan mạnh với các yếu tố như broad_jump_cm, hít đất sit_ups_counts, sit_and_bend_forward_cm, và grip_force, trong khi đó lại tương quan yếu với age, weight_kg, và body_fat. Ngoài ra, cân nặng và tỷ lệ mỡ cơ thể có tương quan dương rất mạnh, còn huyết áp tâm thu và tâm trương cũng liên hệ chặt chẽ với nhau. Điều này nhấn mạnh rằng các yếu tố về sức mạnh, độ bền và dẻo dai là những chỉ số quan trọng đối với hiệu suất, trong khi mỡ cơ thể và tuổi cao thường làm giảm hiệu suất.

• Tóm tắt thống kê

# Tóm tắt thống kê theo giới tính
data %>% group_by(gender) %>%
  summarise(across(c(age, height_cm, weight_kg, body_fat, diastolic, systolic, grip_force, sit_and_bend_forward_cm, sit_ups_counts, broad_jump_cm), mean, na.rm = TRUE))

## Warning: There was 1 warning in `summarise()`.
## ℹ In argument: `across(...)`.
## ℹ In group 1: `gender = "F"`.
## Caused by warning:
## ! The `...` argument of `across()` is deprecated as of dplyr 1.1.0.
## Supply arguments directly to `.fns` through an anonymous function instead.
## 
##   # Previously
##   across(a:b, mean, na.rm = TRUE)
## 
##   # Now
##   across(a:b, \(x) mean(x, na.rm = TRUE))

## # A tibble: 2 × 11
##   gender   age height_cm weight_kg body_fat diastolic systolic grip_force
##   <chr>  <dbl>     <dbl>     <dbl>    <dbl>     <dbl>    <dbl>      <dbl>
## 1 F       37.8      161.      56.9     28.4      75.7     124.       25.9
## 2 M       36.1      173.      73.4     20.2      80.7     134.       43.4
## # ℹ 3 more variables: sit_and_bend_forward_cm <dbl>, sit_ups_counts <dbl>,
## #   broad_jump_cm <dbl>

# Tóm tắt thống kê theo từng nhóm tuổi
library(dplyr)

# Chia age thành 3 nhóm: 20-30, 31-45, 46-64
data <- data %>%
  mutate(age_group = case_when(
    age >= 20 & age <= 34 ~ "20-34",
    age >= 35 & age <= 49 ~ "35-49",
    age >= 50 & age <= 64 ~ "50-64",
    TRUE ~ "Unknown"
  ))
data %>% 
  group_by(age_group) %>%
  summarise(across(c(height_cm, weight_kg, body_fat, diastolic, systolic, grip_force, 
                     sit_and_bend_forward_cm, sit_ups_counts, broad_jump_cm), 
                   mean, na.rm = TRUE))

## # A tibble: 3 × 10
##   age_group height_cm weight_kg body_fat diastolic systolic grip_force
##   <chr>         <dbl>     <dbl>    <dbl>     <dbl>    <dbl>      <dbl>
## 1 20-34          170.      68.1     21.8      77.3     128.       38.4
## 2 35-49          169.      69.0     23.5      80.8     131.       38.0
## 3 50-64          164.      64.2     26.1      80.5     135.       32.7
## # ℹ 3 more variables: sit_and_bend_forward_cm <dbl>, sit_ups_counts <dbl>,
## #   broad_jump_cm <dbl>

# Tóm tắt theo lớp hiệu suất
data %>% group_by(class) %>%
  summarise(across(c(age, height_cm, weight_kg, body_fat, diastolic, systolic, grip_force, sit_and_bend_forward_cm, sit_ups_counts, broad_jump_cm), mean, na.rm = TRUE))

## # A tibble: 4 × 11
##   class   age height_cm weight_kg body_fat diastolic systolic grip_force
##   <chr> <dbl>     <dbl>     <dbl>    <dbl>     <dbl>    <dbl>      <dbl>
## 1 A      35.2      168.      64.4     20.5      77.9     129.       38.6
## 2 B      37.0      169.      66.6     22.0      78.7     131.       38.0
## 3 C      36.6      169.      66.8     22.6      78.6     130.       36.7
## 4 D      38.0      169.      71.8     27.6      80.2     131.       34.8
## # ℹ 3 more variables: sit_and_bend_forward_cm <dbl>, sit_ups_counts <dbl>,
## #   broad_jump_cm <dbl>

• Tóm tắt thống kê theo giới tính: Nam có ưu thế về chiều cao (173.22 cm), cân nặng (73.44 kg), lực nắm tay (43.43 kg), huyết áp (tâm thu 133.81 mmHg, tâm trương 80.67 mmHg), sức bền (44.99 lần gập bụng) và khả năng nhảy xa (211.71 cm). Trong khi đó, nữ vượt trội hơn về tỷ lệ mỡ cơ thể (28.37%), độ dẻo dai (18.83 cm) và có huyết áp thấp hơn (tâm thu 124.03 mmHg, tâm trương 75.68 mmHg). Sự khác biệt này phản ánh đặc điểm sinh lý giữa hai giới.

• Tóm tắt thống kê theo từng nhóm tuổi: Nhóm tuổi 20-34 có ưu thế về chiều cao (170.40 cm), cân nặng (68.07 kg), lực nắm tay (38.43 kg), độ dẻo dai (15.75 cm), sức bền (46.09 lần gập bụng) và khả năng nhảy xa (203.30 cm). Trong khi đó, nhóm tuổi 50-64 có tỷ lệ mỡ cơ thể cao nhất (26.14%), huyết áp cao nhất (tâm thu 134.77 mmHg, tâm trương 80.53 mmHg) nhưng lực nắm tay, sức bền, và khả năng nhảy xa đều giảm. Sự khác biệt này phản ánh ảnh hưởng của tuổi tác đến thể chất và sức khỏe.

• Tóm tắt thống kê theo lớp hiệu suất: Nhóm hiệu suất cao nhất (lớp A) có các chỉ số thể chất vượt trội như lực nắm tay (38.62 kg), khả năng nhảy xa (202.87 cm), sức bền (47.91 lần gập bụng) và độ dẻo dai (21.29 cm). Trong khi đó, nhóm hiệu suất thấp nhất (lớp D) có tỷ lệ mỡ cơ thể cao nhất (27.62%), huyết áp cao nhất (tâm thu 131.10 mmHg, tâm trương 80.17 mmHg), nhưng lực nắm tay, độ dẻo dai, sức bền và khả năng nhảy xa đều thấp hơn đáng kể. Điều này cho thấy hiệu suất giảm có thể liên quan đến tỷ lệ mỡ cơ thể cao và sức khỏe tim mạch kém.

Kiểm định ANOVA cho các biến số học

# Lấy danh sách các biến số học trong dữ liệu
numeric_vars <- names(data)[sapply(data, is.numeric)]  
numeric_vars <- numeric_vars[numeric_vars != "class"] 
numeric_vars <- numeric_vars[numeric_vars != "gender"] 

# Thực hiện kiểm định ANOVA cho từng biến số học
for (var in numeric_vars) {
  formula <- as.formula(paste(var, "~ class"))
  anova_result <- aov(formula, data = data)
  print(paste("ANOVA result for", var))
  print(summary(anova_result))
  cat("\n\n")
}

## [1] "ANOVA result for age"
##                Df  Sum Sq Mean Sq F value   Pr(>F)    
## class           3   13297    4432   24.14 1.43e-15 ***
## Residuals   13121 2409454     184                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for height_cm"
##                Df Sum Sq Mean Sq F value  Pr(>F)    
## class           3   3181  1060.2   15.35 5.8e-10 ***
## Residuals   13121 906538    69.1                    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for weight_kg"
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## class           3   95475   31825   245.4 <2e-16 ***
## Residuals   13121 1701275     130                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for body_fat"
##                Df Sum Sq Mean Sq F value Pr(>F)    
## class           3  93643   31214   716.5 <2e-16 ***
## Residuals   13121 571649      44                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for diastolic"
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## class           3    8788  2929.5   26.95 <2e-16 ***
## Residuals   13121 1426084   108.7                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for systolic"
##                Df  Sum Sq Mean Sq F value   Pr(>F)    
## class           3    6269  2089.6   10.12 1.17e-06 ***
## Residuals   13121 2708611   206.4                     
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for grip_force"
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## class           3   27981    9327   85.32 <2e-16 ***
## Residuals   13121 1434366     109                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for sit_and_bend_forward_cm"
##                Df Sum Sq Mean Sq F value Pr(>F)    
## class           3 319289  106430    2729 <2e-16 ***
## Residuals   13121 511694      39                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for sit_ups_counts"
##                Df  Sum Sq Mean Sq F value Pr(>F)    
## class           3  556290  185430    1167 <2e-16 ***
## Residuals   13121 2085660     159                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## [1] "ANOVA result for broad_jump_cm"
##                Df   Sum Sq Mean Sq F value Pr(>F)    
## class           3  1415679  471893   330.9 <2e-16 ***
## Residuals   13121 18710769    1426                   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

p-value của các biến định lượng đều nhỏ hơn 0.05. Do đó, với mỗi biến định lượng thì có ít nhất một nhóm hiệu suất khác với trung bình các nhóm còn lại.

Xây dựng và đánh giá mô hình hồi quy logistic đa biến

Multinominal

• Kiểm tra dữ liệu

# Kiểm tra giá trị thiếu
colSums(is.na(data))

##                     age                  gender               height_cm 
##                       0                       0                       0 
##               weight_kg                body_fat               diastolic 
##                       0                       0                       0 
##                systolic              grip_force sit_and_bend_forward_cm 
##                       0                       0                       0 
##          sit_ups_counts           broad_jump_cm                   class 
##                       0                       0                       0 
##               age_group 
##                       0

# Chuyển đổi các biến phân loại thành factor
data$gender <- as.factor(data$gender)
data$class <- as.factor(data$class)

# Kiểm tra lại cấu trúc sau khi chuyển đổi
str(data)

## gropd_df [13,125 × 13] (S3: grouped_df/tbl_df/tbl/data.frame)
##  $ age                    : num [1:13125] 27 25 31 32 28 36 42 33 54 28 ...
##  $ gender                 : Factor w/ 2 levels "F","M": 2 2 2 2 2 1 1 2 2 2 ...
##  $ height_cm              : num [1:13125] 172 165 180 174 174 ...
##  $ weight_kg              : num [1:13125] 75.2 55.8 78 71.1 67.7 ...
##  $ body_fat               : num [1:13125] 21.3 15.7 20.1 18.4 17.1 22 32.2 36.9 27.6 14.4 ...
##  $ diastolic              : num [1:13125] 80 77 92 76 70 64 72 84 85 81 ...
##  $ systolic               : num [1:13125] 130 126 152 147 127 119 135 137 165 156 ...
##  $ grip_force             : num [1:13125] 54.9 36.4 44.8 41.4 43.5 23.8 22.7 45.9 40.4 57.9 ...
##  $ sit_and_bend_forward_cm: num [1:13125] 18.4 16.3 12 15.2 27.1 21 0.8 12.3 18.6 12.1 ...
##  $ sit_ups_counts         : num [1:13125] 60 53 49 53 45 27 18 42 34 55 ...
##  $ broad_jump_cm          : num [1:13125] 217 229 181 219 217 153 146 234 148 213 ...
##  $ class                  : Factor w/ 4 levels "A","B","C","D": 3 1 3 2 2 2 4 2 3 2 ...
##  $ age_group              : chr [1:13125] "20-34" "20-34" "20-34" "20-34" ...
##  - attr(*, "groups")= tibble [4 × 2] (S3: tbl_df/tbl/data.frame)
##   ..$ class: Factor w/ 4 levels "A","B","C","D": 1 2 3 4
##   ..$ .rows: list<int> [1:4] 
##   .. ..$ : int [1:3303] 2 11 16 18 22 30 33 38 39 45 ...
##   .. ..$ : int [1:3285] 4 5 6 8 10 19 20 24 35 47 ...
##   .. ..$ : int [1:3283] 1 3 9 13 14 15 17 21 25 26 ...
##   .. ..$ : int [1:3254] 7 12 23 27 28 34 37 40 41 42 ...
##   .. ..@ ptype: int(0) 
##   ..- attr(*, ".drop")= logi TRUE

• Biểu đồ phân phối biến phân loại class

library(plotly)

## 
## Attaching package: 'plotly'

## The following object is masked from 'package:ggplot2':
## 
##     last_plot

## The following object is masked from 'package:stats':
## 
##     filter

## The following object is masked from 'package:graphics':
## 
##     layout

bieu_do_class <- ggplot(data, aes(x = class)) +
  geom_bar(aes(fill = class), color = "black", stat = "count") +
  labs(title = "Bar Plot of Class Counts", x = "Class", y = "Frequency") +
  theme_minimal()


ggplotly(bieu_do_class)

• Tách dữ liệu thành tập huấn luyện và kiểm tra

set.seed(123)
train_index <- sample(1:nrow(data), 0.7 * nrow(data))
train_data <- data[train_index, ]
test_data <- data[-train_index, ]

• Xây dựng mô hình multinominal logistic

library(nnet)  # Thư viện chứa hàm multinom()

data_md <- multinom(formula = class ~ age + weight_kg + height_cm + body_fat + diastolic + systolic + grip_force + sit_and_bend_forward_cm + sit_ups_counts + broad_jump_cm + gender,
data = train_data, maxit = 1500)

## # weights:  52 (36 variable)
## initial  value 12735.886296 
## iter  10 value 9573.635266
## iter  20 value 9496.969819
## iter  30 value 9367.270822
## iter  40 value 7818.617886
## iter  50 value 7816.882402
## iter  50 value 7816.882349
## iter  50 value 7816.882348
## final  value 7816.882348 
## converged

• Dự đoán và đánh giá trên tập kiểm tra

# Dự đoán trên tập kiểm tra
pred <- predict(data_md, newdata = test_data)

# Độ chính xác
accuracy <- mean(pred == test_data$class)
print(paste("Accuracy:", accuracy))

## [1] "Accuracy: 0.619857795835449"

# Ma trận nhầm lẫn
conf_matrix <- table(Predicted = pred, Actual = test_data$class)
print(conf_matrix)

##          Actual
## Predicted   A   B   C   D
##         A 739 205  72  12
##         B 259 464 192  60
##         C  22 289 506 150
##         D   0  39 197 732

• Kiểm tra giá trị p-value

z_values <- summary(data_md)$coefficients / summary(data_md)$standard.errors
p_values <- (1 - pnorm(abs(z_values), 0, 1)) * 2
print(p_values)

##   (Intercept) age weight_kg    height_cm     body_fat    diastolic   systolic
## B           0   0         0 1.366443e-01 8.888073e-01 1.831951e-01 0.27833114
## C           0   0         0 5.829557e-03 6.591629e-01 5.254671e-03 0.04935171
## D           0   0         0 2.781775e-12 4.778272e-09 7.330586e-05 0.02393082
##   grip_force sit_and_bend_forward_cm sit_ups_counts broad_jump_cm genderM
## B          0                       0              0  6.661338e-16       0
## C          0                       0              0  0.000000e+00       0
## D          0                       0              0  0.000000e+00       0

Dùng nhóm A làm nhóm tham chiếu cho các nhóm B, C và D so sánh, đồng thời nhìn vào bảng số liệu ta có thể thấy rằng:

• Biến height_cm trong nhóm B có p_val > 0.05.

• Biến body_fat trong nhóm B, C đều có p_val > 0.05

• Biến diastolic trong nhóm B, C có p_val > 0.05.

• Biến systolic trong nhóm B đều có p_val > 0.05.

Qua việc xử lý như trên, có thể thấy biến body_fat và biến diastolic không có ảnh hưởng gì nhiều đến hiệu suất của mô hình phân loại biến class.

Còn các biến còn lại cũng không có ảnh hưởng hoặc sẽ ảnh hưởng ít đến quá trình thống kê và đưa ra dự đoán nên có thể loại bỏ các biến này trong mô hình.

Nhìn chung, các biến này đều liên quan đến nhóm B.

• Huấn luyện lại mô hình multinomial logistic regression

library(nnet)
data_md <- multinom(
  formula = class ~ age + weight_kg + grip_force + 
    sit_and_bend_forward_cm + sit_ups_counts + broad_jump_cm + gender,
  data = train_data,
  maxit = 1500
)

## # weights:  36 (24 variable)
## initial  value 12735.886296 
## iter  10 value 10000.438877
## iter  20 value 9897.648406
## iter  30 value 7904.115215
## final  value 7903.783114 
## converged

# Dự đoán trên tập kiểm tra
pred <- predict(data_md, newdata = test_data)

# Tính độ chính xác
accuracy <- mean(pred == test_data$class)
print(paste("Accuracy:", accuracy))

## [1] "Accuracy: 0.618080243778568"

# Tạo ma trận nhầm lẫn
conf_matrix <- table(Predicted = pred, Actual = test_data$class)

# Kiểm tra xem conf_matrix có đúng định dạng không
if (!is.matrix(conf_matrix)) {
  conf_matrix <- as.matrix(conf_matrix)
}

# In ma trận nhầm lẫn
print("Confusion Matrix:")

## [1] "Confusion Matrix:"

print(conf_matrix)

##          Actual
## Predicted   A   B   C   D
##         A 739 202  75  12
##         B 261 469 196  55
##         C  20 291 507 168
##         D   0  35 189 719

Sau khi loại bỏ các biến không ảnh hưởng thì hiệu suất là 0.618, so với hiệu suất lúc ban đầu là 0.619, hiệu suất của mô hình không giảm nhiều so với ban đầu. Điều này cho thấy các biến đã bị loại bỏ thật sự không ảnh hưởng đáng kể tới hiệu suất của mô hình.

• Hàm đánh giá đa lớp

# Hàm đánh giá
eval_multi_class <- function(conf_matrix) {
  cc <- sum(diag(conf_matrix)) # Số dự đoán đúng
  sc <- sum(conf_matrix)       # Tổng số mẫu
  
  pp <- colSums(conf_matrix) # Tổng dự đoán cho mỗi lớp
  tt <- rowSums(conf_matrix) # Tổng thực tế cho mỗi lớp
  
  precision <- diag(conf_matrix) / pp
  recall <- diag(conf_matrix) / tt
  f1_scores <- 2 * precision * recall / (precision + recall)
  
  macro_precision <- mean(precision, na.rm = TRUE)
  macro_recall <- mean(recall, na.rm = TRUE)
  macro_f1 <- mean(f1_scores, na.rm = TRUE)
  
  accuracy <- cc / sc
  
  expected <- (pp * tt) / sc
  kappa <- (cc - sum(expected)) / (sc - sum(expected))
  
  return(list(
    Precision = precision,
    Recall = recall,
    F1_Scores = f1_scores,
    Macro_Precision = macro_precision,
    Macro_Recall = macro_recall,
    Macro_F1 = macro_f1,
    Accuracy = accuracy,
    Kappa = kappa
  ))
}

# Gọi hàm đánh giá
results <- eval_multi_class(conf_matrix)

# Hiển thị kết quả đánh giá
print("Evaluation Metrics:")

## [1] "Evaluation Metrics:"

print(results)

## $Precision
##         A         B         C         D 
## 0.7245098 0.4704112 0.5243020 0.7536688 
## 
## $Recall
##         A         B         C         D 
## 0.7188716 0.4780836 0.5141988 0.7624602 
## 
## $F1_Scores
##         A         B         C         D 
## 0.7216797 0.4742164 0.5192012 0.7580390 
## 
## $Macro_Precision
## [1] 0.6182229
## 
## $Macro_Recall
## [1] 0.6184035
## 
## $Macro_F1
## [1] 0.6182841
## 
## $Accuracy
## [1] 0.6180802
## 
## $Kappa
## [1] 0.4906537

Mô hình hoạt động tốt với các lớp A và D, với F1-Score lần lượt là 0.721 và 0.758, nhưng hiệu suất còn hạn chế ở các lớp B (F1-Score 0.474) và C (F1-Score 0.519). Độ chính xác tổng thể đạt 61.80%, các chỉ số tổng hợp như Macro Precision (0.618), Macro Recall (0.618), và Macro F1 (0.618) cho thấy mô hình duy trì sự cân đối nhưng chưa thực sự vượt trội. Chỉ số Kappa (0.490) phản ánh mức độ đồng thuận trung bình giữa mô hình và dữ liệu thực tế, tốt hơn dự đoán ngẫu nhiên nhưng vẫn cần cải thiện.

• Vẽ biểu đồ tương quan cho ma trận nhầm lẫn

library(ggplot2)
library(reshape2)

# Heatmap ma trận nhầm lẫn
conf_matrix <- melt(as.matrix(conf_matrix))
ggplot(conf_matrix, aes(x = Predicted, y = Actual, fill = value)) +
  geom_tile() +
  geom_text(aes(label = value), color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  labs(title = "Confusion Matrix Heatmap", x = "Predicted", y = "Actual")

Heatmap của ma trận nhầm lẫn cho thấy mô hình hoạt động tốt nhất trên Class A và Class D, với số lượng phân loại đúng lần lượt là 739 và 719. Tuy nhiên, vẫn xảy ra nhầm lẫn đáng kể, đặc biệt giữa Class A với Class B (202 mẫu nhầm) và giữa Class C với Class B (291 mẫu nhầm). Class B và Class C có sự nhầm lẫn lẫn nhau khá nhiều, cho thấy các đặc trưng phân biệt giữa hai lớp này chưa rõ ràng. Điều này có thể là nguyên nhân dẫn đến sự chênh lệch trong hiệu suất phân loại giữa các lớp.

LDA

• Huấn luyện mô hình LDA

library(MASS)

## 
## Attaching package: 'MASS'

## The following object is masked from 'package:plotly':
## 
##     select

## The following object is masked from 'package:dplyr':
## 
##     select

data_md <- lda(formula = class ~ age + weight_kg + height_cm + body_fat +systolic + diastolic + grip_force + sit_and_bend_forward_cm + sit_ups_counts + broad_jump_cm + gender,
data = train_data, maxit = 1500)

data_md

## Call:
## lda(class ~ age + weight_kg + height_cm + body_fat + systolic + 
##     diastolic + grip_force + sit_and_bend_forward_cm + sit_ups_counts + 
##     broad_jump_cm + gender, data = train_data, maxit = 1500)
## 
## Prior probabilities of groups:
##         A         B         C         D 
## 0.2485033 0.2490476 0.2520954 0.2503538 
## 
## Group means:
##        age weight_kg height_cm body_fat systolic diastolic grip_force
## A 35.14455  64.44451  167.9035 20.45569  129.463  77.96163   38.69708
## B 36.89729  66.52468  168.6502 21.92157  130.281  78.44580   38.01610
## C 36.52763  66.74383  169.2620 22.47904  129.674  78.47064   36.66625
## D 37.98783  71.69750  168.6080 27.64033  130.830  80.03883   34.69743
##   sit_and_bend_forward_cm sit_ups_counts broad_jump_cm   genderM
## A               21.303316       48.14717      203.2501 0.5606658
## B               17.382740       42.73837      195.6967 0.6464161
## C               14.323402       38.86010      189.4685 0.6727116
## D                7.924117       30.01478      174.1665 0.6565217
## 
## Coefficients of linear discriminants:
##                                  LD1          LD2          LD3
## age                     -0.046551687  0.015719167  0.034160234
## weight_kg                0.046848417  0.092133536 -0.003805877
## height_cm               -0.007376237 -0.110431258 -0.026165844
## body_fat                 0.025086802  0.032983744  0.087662824
## systolic                -0.003026230 -0.001282414  0.014818755
## diastolic                0.006385492  0.004931975 -0.035825151
## grip_force              -0.044335238  0.030505319  0.016434001
## sit_and_bend_forward_cm -0.103802771 -0.031192913  0.059784637
## sit_ups_counts          -0.071621443  0.025692202  0.001128416
## broad_jump_cm           -0.006388751  0.016436049 -0.002191084
## genderM                  1.297979683 -2.540723971  2.422138229
## 
## Proportion of trace:
##    LD1    LD2    LD3 
## 0.9799 0.0186 0.0015

• Dự đoán và đánh giá trên tập kiểm tra

# Dự đoán trên dữ liệu
predictions <- predict(data_md, newdata = test_data)

# Tạo ma trận nhầm lẫn
conf_matrix <- table(Predicted = predictions$class, Actual = test_data$class)
print(conf_matrix)

##          Actual
## Predicted   A   B   C   D
##         A 725 217  77  13
##         B 265 440 189  67
##         C  30 316 546 174
##         D   0  24 155 700

# Tính toán độ chính xác
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
cat("Accuracy:", accuracy, "\n")

## Accuracy: 0.6122397

• Hệ số của mô hình LDA

data_md$scaling  # Các hệ số của từng biến độc lập trong các hàm phân biệt

##                                  LD1          LD2          LD3
## age                     -0.046551687  0.015719167  0.034160234
## weight_kg                0.046848417  0.092133536 -0.003805877
## height_cm               -0.007376237 -0.110431258 -0.026165844
## body_fat                 0.025086802  0.032983744  0.087662824
## systolic                -0.003026230 -0.001282414  0.014818755
## diastolic                0.006385492  0.004931975 -0.035825151
## grip_force              -0.044335238  0.030505319  0.016434001
## sit_and_bend_forward_cm -0.103802771 -0.031192913  0.059784637
## sit_ups_counts          -0.071621443  0.025692202  0.001128416
## broad_jump_cm           -0.006388751  0.016436049 -0.002191084
## genderM                  1.297979683 -2.540723971  2.422138229

data_md$svd

## [1] 72.043808  9.918424  2.838856

• Hàm đánh giá đa lớp

# Gọi hàm đánh giá
results <- eval_multi_class(conf_matrix)

# Hiển thị kết quả đánh giá
print("Evaluation Metrics:")

## [1] "Evaluation Metrics:"

print(results)

## $Precision
##         A         B         C         D 
## 0.7107843 0.4413240 0.5646329 0.7337526 
## 
## $Recall
##         A         B         C         D 
## 0.7025194 0.4578564 0.5121951 0.7963595 
## 
## $F1_Scores
##         A         B         C         D 
## 0.7066277 0.4494382 0.5371372 0.7637752 
## 
## $Macro_Precision
## [1] 0.6126234
## 
## $Macro_Recall
## [1] 0.6172326
## 
## $Macro_F1
## [1] 0.6142446
## 
## $Accuracy
## [1] 0.6122397
## 
## $Kappa
## [1] 0.4828447

Kết quả đánh giá mô hình LDA cho thấy độ chính xác tổng thể đạt 61.22%, phản ánh hiệu suất trung bình trong việc phân loại dữ liệu. Macro Precision, Recall, và F1-Score lần lượt đạt 61.26%, 61.72%, và 61.42%, cho thấy sự nhất quán giữa các chỉ số. Tuy nhiên, hiệu suất phân loại khác biệt đáng kể giữa các lớp. Lớp A và D có F1-Score lần lượt là 70.66% và 76.37%, thể hiện khả năng phân loại tốt. Trong khi đó, lớp B có hiệu suất thấp nhất với F1-Score chỉ đạt 44.94%, cho thấy mô hình gặp khó khăn trong việc phân biệt lớp này. Giá trị Kappa là 0.48, chỉ ra rằng mô hình có mức độ phù hợp trung bình so với một mô hình dự đoán ngẫu nhiên.

• Vẽ biểu đồ tương quan cho ma trận nhầm lẫn

library(ggplot2)
library(reshape2)

# Heatmap ma trận nhầm lẫn
conf_matrix <- melt(as.matrix(conf_matrix))
ggplot(conf_matrix, aes(x = Predicted, y = Actual, fill = value)) +
  geom_tile() +
  geom_text(aes(label = value), color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  labs(title = "Confusion Matrix Heatmap", x = "Predicted", y = "Actual")

Heatmap của ma trận nhầm lẫn cho thấy mô hình hoạt động tốt nhất trên Class A và Class D, với số lượng phân loại đúng lần lượt là 725 và 700. Tuy nhiên, vẫn xảy ra nhầm lẫn đáng kể, đặc biệt giữa Class A với Class B (265 mẫu nhầm) và giữa Class C với Class B (189 mẫu nhầm). Class B và Class C có sự nhầm lẫn lẫn nhau khá nhiều, cho thấy các đặc trưng phân biệt giữa hai lớp này chưa rõ ràng. Điều này có thể là nguyên nhân dẫn đến sự chênh lệch trong hiệu suất phân loại giữa các lớp.

• Kiểm tra tính đồng nhất của ma trận hiệp phương sai (Box’s M Test)

x <- train_data[, c("age", "weight_kg", "height_cm", "body_fat", "grip_force", "diastolic", "systolic",
                    "sit_and_bend_forward_cm", "sit_ups_counts", 
                    "broad_jump_cm")]
x$gender <- as.numeric(factor(train_data$gender))  # Chuyển gender thành số (1 cho M, 0 cho F)

# Biến phụ thuộc (nhóm phân loại)
y <- train_data$class

# Thực hiện Box's M Test
library(biotools)

## ---
## biotools version 4.2

box_m_test <- boxM(as.matrix(x), grouping = y)

# In kết quả
print(box_m_test)

## 
##  Box's M-test for Homogeneity of Covariance Matrices
## 
## data:  as.matrix(x)
## Chi-Sq (approx.) = 5390.6, df = 198, p-value < 2.2e-16

Do p_vals < 0.05 nên không thể cho rằng giả định các nhóm đồng nhất phương sai với nhau. Việc không thể chắc chắn rằng các nhóm đồng nhất phương sai với nhau thì có thể gây ra ảnh hưởng về hiệu quả dự đoán của mô hình LDA.

• Phân tích phân phối các biến độc lập

library(ggplot2)
library(reshape2)

# Danh sách các biến độc lập
numeric_vars <- c("age", "height_cm", "weight_kg", "body_fat", "diastolic", 
                  "systolic", "grip_force", "sit_and_bend_forward_cm", 
                  "sit_ups_counts", "broad_jump_cm")

# Chuyển đổi dữ liệu từ wide format sang long format để dễ dàng vẽ
data_long <- melt(data, id.vars = "class", measure.vars = numeric_vars)

# Vẽ phân phối cho tất cả các biến độc lập
ggplot(data_long, aes(x = value, fill = as.factor(class))) + 
  geom_histogram(position = "dodge", bins = 30) + 
  facet_wrap(~ variable, scales = "free") + # Tạo các biểu đồ con cho mỗi biến
  labs(title = "Distribution of Independent Variables by Class", x = "Value", y = "Count") +
  theme_minimal()

Qua việc kiểm tra phân phối chuẩn của từng nhóm dữ liệu loại A, B, C và D cho thấy rằng biến age, grip_force và broad_jump_cm không tuân theo phân phối chuẩn nên cũng ảnh hưởng đến việc phân loại của mô hình.

Do việc sử dụng LDA có thể không đảm bảo độ chính xác, nên việc chuyển sang QDA là hợp lý hơn, vì QDA không yêu cầu giả định về sự đồng nhất phương sai và phân phối của dữ liệu.

QDA

• Huấn luyện mô hình QDA

qda_model <- qda(class ~ age + weight_kg + height_cm + body_fat + diastolic + systolic + grip_force + sit_and_bend_forward_cm + 
                 sit_ups_counts + broad_jump_cm + gender, data = train_data)

• Dự đoán trên tập kiểm tra

# Dự đoán trên tập kiểm tra
pred_qda <- predict(qda_model, newdata = test_data)$class
# Tạo ma trận nhầm lẫn
conf_matrix_qda <- table(Predicted = pred_qda, Actual = test_data$class)
print(conf_matrix_qda)

##          Actual
## Predicted   A   B   C   D
##         A 773 185  65  10
##         B 234 564 206  61
##         C   9 223 620 160
##         D   4  25  76 723

# Đánh giá mô hình
accuracy_qda <- mean(pred_qda == test_data$class)
cat("Accuracy:", accuracy_qda)

## Accuracy: 0.6805485

qda_model

## Call:
## qda(class ~ age + weight_kg + height_cm + body_fat + diastolic + 
##     systolic + grip_force + sit_and_bend_forward_cm + sit_ups_counts + 
##     broad_jump_cm + gender, data = train_data)
## 
## Prior probabilities of groups:
##         A         B         C         D 
## 0.2485033 0.2490476 0.2520954 0.2503538 
## 
## Group means:
##        age weight_kg height_cm body_fat diastolic systolic grip_force
## A 35.14455  64.44451  167.9035 20.45569  77.96163  129.463   38.69708
## B 36.89729  66.52468  168.6502 21.92157  78.44580  130.281   38.01610
## C 36.52763  66.74383  169.2620 22.47904  78.47064  129.674   36.66625
## D 37.98783  71.69750  168.6080 27.64033  80.03883  130.830   34.69743
##   sit_and_bend_forward_cm sit_ups_counts broad_jump_cm   genderM
## A               21.303316       48.14717      203.2501 0.5606658
## B               17.382740       42.73837      195.6967 0.6464161
## C               14.323402       38.86010      189.4685 0.6727116
## D                7.924117       30.01478      174.1665 0.6565217

Trong 4 nhóm, các biến age, height_cm, systolic, diastolic, và gender không có sự chênh lệch đáng kể (không quá 5%). Vì vậy, các biến này được loại bỏ để thực hiện lại mô hình.

• Chạy lại mô hình sau khi bỏ đi 5 biến nêu trên

qda_model <- qda(class ~ weight_kg + body_fat + grip_force + sit_and_bend_forward_cm + 
                 sit_ups_counts + broad_jump_cm, data = train_data)
# Dự đoán trên tập kiểm tra
pred_qda <- predict(qda_model, newdata = test_data)$class

# Tạo ma trận nhầm lẫn
conf_matrix_qda <- table(Predicted = pred_qda, Actual = test_data$class)
print(conf_matrix_qda)

##          Actual
## Predicted   A   B   C   D
##         A 747 254 102  17
##         B 223 410 208  59
##         C  44 287 539 176
##         D   6  46 118 702

# Đánh giá mô hình
accuracy_qda <- mean(pred_qda == test_data$class)
cat("Accuracy:", accuracy_qda)

## Accuracy: 0.6089385

Kết quả chạy lại mô hình cho thấy việc loại bỏ các biến trên không hoàn toàn chính xác:

• Biến age: Mặc dù không có sự chênh lệch trong 4 nhóm A, B, C và D và có tương quan thấp với biến class, nhưng đồ thị phân phối cho thấy biến này không tuân theo phân phối chuẩn. Hơn nữa, trong đồ thị tương quan, age có mối quan hệ mạnh mẽ với hai biến sit_ups_counts và broad_jump_cm. Do đó, biến age đóng vai trò bổ trợ và không nên loại bỏ.

• Biến height_cm: Phân phối chuẩn nhưng có ảnh hưởng mạnh đến các biến khác trong ma trận tương quan và ít ảnh hưởng đến biến class. Điều này cho thấy khả năng xảy ra đa cộng tuyến, nên có thể loại bỏ biến này.

• Biến systolic và diastolic: Cả hai không có tương quan đáng kể với các biến khác và không có sự chênh lệch trung bình trong 4 nhóm A, B, C và D. Vì vậy, việc loại bỏ hai biến này là hợp lý.

• Biến gender: Là biến định tính phân biệt nam và nữ, cần giữ lại do biểu đồ heatmap cho thấy gender có ảnh hưởng đến các biến khác, tương tự như vai trò của biến age.

• Thử nghiệm các biến độc lập khác với QDA

qda_model <- qda(class ~ age + weight_kg + body_fat + grip_force + sit_and_bend_forward_cm + 
                 sit_ups_counts + broad_jump_cm + gender, data = train_data)
# Dự đoán trên tập kiểm tra
pred_qda <- predict(qda_model, newdata = test_data)$class

# Tạo ma trận nhầm lẫn
conf_matrix_qda <- table(Predicted = pred_qda, Actual = test_data$class)
print(conf_matrix_qda)

##          Actual
## Predicted   A   B   C   D
##         A 768 184  63  12
##         B 241 565 217  58
##         C   7 225 606 159
##         D   4  23  81 725

# Đánh giá mô hình
accuracy_qda <- mean(pred_qda == test_data$class)
cat("Accuracy:", accuracy_qda)

## Accuracy: 0.6764855

Kết quả cho thấy khi sử dụng tất cả các biến, độ chính xác đạt 0,680. Tuy nhiên, sau khi loại bỏ ba biến height_cm, diastolic, và systolic, độ chính xác giảm nhẹ xuống 0,676. Điều này chứng tỏ các biến này không ảnh hưởng nhiều đến hiệu suất của mô hình.

Ngoài ra, khi xét biến body_fat trong mô hình Multinomial Logistic, kết quả cho thấy biến này không có vai trò quan trọng trong việc phân loại các giá trị của biến class. Do đó, việc loại bỏ biến body_fat cũng có thể được cân nhắc để đơn giản hóa mô hình.

• Chạy lại mô hình sau khi loại bỏ biến body_fat

qda_model <- qda(class ~ age + weight_kg + grip_force + sit_and_bend_forward_cm + 
                 sit_ups_counts + broad_jump_cm + gender, data = train_data)
# Dự đoán trên tập kiểm tra
pred_qda <- predict(qda_model, newdata = test_data)$class

# Tạo ma trận nhầm lẫn
conf_matrix_qda <- table(Predicted = pred_qda, Actual = test_data$class)
print(conf_matrix_qda)

##          Actual
## Predicted   A   B   C   D
##         A 772 183  64   9
##         B 239 559 211  47
##         C   6 232 602 181
##         D   3  23  90 717

# Đánh giá mô hình
accuracy_qda <- mean(pred_qda == test_data$class)
cat("Accuracy:", accuracy_qda)

## Accuracy: 0.6729304

Hiệu suất của mô hình chỉ giảm 0,4%, điều này cho thấy biến body_fat có thể gặp vấn đề đa cộng tuyến giống như biến height_cm và thực sự không ảnh hưởng đáng kể đến mô hình. Vì vậy, việc loại bỏ biến body_fat là hợp lý.

• Hàm đánh giá đa lớp

# Gọi hàm đánh giá
results <- eval_multi_class(conf_matrix_qda)

# Hiển thị kết quả đánh giá
print("Evaluation Metrics:")

## [1] "Evaluation Metrics:"

print(results)

## $Precision
##         A         B         C         D 
## 0.7568627 0.5606820 0.6225440 0.7515723 
## 
## $Recall
##         A         B         C         D 
## 0.7509728 0.5293561 0.5896180 0.8607443 
## 
## $F1_Scores
##         A         B         C         D 
## 0.7539062 0.5445689 0.6056338 0.8024622 
## 
## $Macro_Precision
## [1] 0.6729153
## 
## $Macro_Recall
## [1] 0.6826728
## 
## $Macro_F1
## [1] 0.6766428
## 
## $Accuracy
## [1] 0.6729304
## 
## $Kappa
## [1] 0.5636664

Độ chính xác tổng thể của mô hình đạt 67.29%, phản ánh hiệu suất khá tốt trong việc phân loại dữ liệu. Macro Precision, Recall, và F1-Score lần lượt đạt 67.29%, 68.27%, và 67.66%, cho thấy sự nhất quán giữa các chỉ số tổng thể. Xét từng lớp, lớp A và D có F1-Score cao (lần lượt là 75.39% và 80.25%), thể hiện khả năng phân loại tốt. Trong khi đó, lớp B có F1-Score thấp nhất (54.46%), cho thấy mô hình gặp khó khăn trong việc phân biệt lớp này. Lớp C có F1-Score trung bình (60.56%), phản ánh hiệu suất phân loại ổn định hơn. Giá trị Kappa đạt 0.56, chỉ ra rằng mô hình có mức độ phù hợp khá cao so với một mô hình dự đoán ngẫu nhiên.

• Vẽ biểu đồ tương quan cho ma trận nhầm lẫn

library(ggplot2)
library(reshape2)

# Heatmap ma trận nhầm lẫn
conf_matrix <- melt(as.matrix(conf_matrix_qda))
ggplot(conf_matrix, aes(x = Predicted, y = Actual, fill = value)) +
  geom_tile() +
  geom_text(aes(label = value), color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  labs(title = "Confusion Matrix Heatmap", x = "Predicted", y = "Actual")

Heatmap của ma trận nhầm lẫn cho thấy mô hình hoạt động tốt nhất trên Class A và Class D, với số lượng phân loại đúng lần lượt là 772 và 717. Tuy nhiên, vẫn xảy ra nhầm lẫn đáng kể, đặc biệt giữa Class A với Class B (239 mẫu nhầm) và giữa Class B với Class C (232 mẫu nhầm). Class B và Class C có sự nhầm lẫn lẫn nhau khá nhiều, cho thấy các đặc trưng phân biệt giữa hai lớp này chưa rõ ràng. Điều này có thể là nguyên nhân dẫn đến sự chênh lệch trong hiệu suất phân loại giữa các lớp.

RandomForest

Qua việc đọc slide và tìm hiểu tài liệu, nhóm em thấy rằng phương pháp phân nhóm khi sử dụng model “RandomForest” phù hợp với bài toán (có biến định tính lẫn biến định lương, nhiều biến đầu vào,…).

• Huấn luyện mô hình Random Forest

library(randomForest)

## randomForest 4.7-1.2

## Type rfNews() to see new features/changes/bug fixes.

## 
## Attaching package: 'randomForest'

## The following object is masked from 'package:dplyr':
## 
##     combine

## The following object is masked from 'package:ggplot2':
## 
##     margin

# Xây dựng mô hình Random Forest với tất cả các biến
rf_model <- randomForest(class ~ ., data = train_data, ntree = 500, mtry = 3, importance = TRUE)

# In kết quả mô hình
print(rf_model)

## 
## Call:
##  randomForest(formula = class ~ ., data = train_data, ntree = 500,      mtry = 3, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 26.27%
## Confusion matrix:
##      A    B    C    D class.error
## A 1970  282   26    5   0.1371003
## B  525 1371  315   77   0.4007867
## C  199  457 1539  121   0.3354922
## D   27  120  259 1894   0.1765217

• Dự đoán và đánh giá trên tập kiểm tra

# Dự đoán
rf_predictions <- predict(rf_model, newdata = test_data)

# Đánh giá mô hình
conf_matrix <- table(Predicted = rf_predictions, Actual = test_data$class)
print(conf_matrix)

##          Actual
## Predicted   A   B   C   D
##         A 859 202  81  11
##         B 150 634 176  50
##         C   8 120 643 108
##         D   3  41  67 785

# Tính độ chính xác
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(paste("Accuracy:", accuracy))

## [1] "Accuracy: 0.741747079735907"

• Hàm đánh giá đa lớp

# Gọi hàm đánh giá
results <- eval_multi_class(conf_matrix)

# Hiển thị kết quả đánh giá
print("Evaluation Metrics:")

## [1] "Evaluation Metrics:"

print(results)

## $Precision
##         A         B         C         D 
## 0.8421569 0.6359077 0.6649431 0.8228512 
## 
## $Recall
##         A         B         C         D 
## 0.7450130 0.6277228 0.7315131 0.8761161 
## 
## $F1_Scores
##         A         B         C         D 
## 0.7906121 0.6317887 0.6966414 0.8486486 
## 
## $Macro_Precision
## [1] 0.7414647
## 
## $Macro_Recall
## [1] 0.7450912
## 
## $Macro_F1
## [1] 0.7419227
## 
## $Accuracy
## [1] 0.7417471
## 
## $Kappa
## [1] 0.6553414

Mô hình có khả năng dự đoán chính xác 74.17% ở mức khá và có thể chấp nhận được. Nhìn vào Precision Recall, F1 và Accuracy các chỉ số khá là tương đương nhau chứng tỏ hiệu suất mô hình khá ổn định. Giá trị Kappa là 0.655 tuy không cao nhưng cũng có thể chấp nhận được đối với mô hình khi phân loại trong bài toán trên. Phân tích từng lớp cho thấy khả năng phân loại hiệu suất A và D khá tốt nhưng ở lớp B và C thì giá trị khá thấp ( có thể do nhiều nguyên nhân như dữ liệu dành cho lớp B và C chưa đủ tốt).

• Kiểm tra sự quan trọng của các biến đầu vào

importance(rf_model)

##                                 A          B          C          D
## age                      62.55917  54.253264  48.903472   8.598975
## gender                   35.49852  28.909149  24.258033  11.838716
## height_cm                30.63456  21.570471  19.577899  18.226819
## weight_kg                58.33983  31.281657  45.188508  39.997104
## body_fat                 42.73943  26.085809  63.431190  53.168123
## diastolic                10.43764   2.185702   6.348282   6.037018
## systolic                 12.61870   5.995448   8.877820   5.501251
## grip_force               68.11862  42.481408  23.395533  11.067115
## sit_and_bend_forward_cm 248.76446 104.620556 105.636503 125.967213
## sit_ups_counts          164.96260  84.765636  67.347935  67.733741
## broad_jump_cm            68.01216  27.198388  16.653955  20.608215
## age_group                32.71239  31.594038  27.427462   3.560962
##                         MeanDecreaseAccuracy MeanDecreaseGini
## age                                 73.15064        507.00245
## gender                              39.07735         96.44737
## height_cm                           43.02584        435.54180
## weight_kg                           77.52990        609.18772
## body_fat                            74.88378        678.29054
## diastolic                           12.59468        337.83095
## systolic                            17.17488        354.04253
## grip_force                          74.54972        553.58873
## sit_and_bend_forward_cm            200.76782       1718.41726
## sit_ups_counts                     131.09241        958.16350
## broad_jump_cm                       65.01324        519.23668
## age_group                           37.55080        120.53625

Dựa trên kết quả kiểm tra mức độ quan trọng của các biến, các biến “diastolic”, “systolic”, “height_cm” và “gender” đều không có tác động đáng kể hoặc ảnh hưởng rất ít đến hiệu suất của mô hình Random Forest và các mô hình khác trước đó. Do đó, quyết định được đưa ra là loại bỏ 4 biến này nhằm giảm sự phức tạp của mô hình, qua đó cải thiện hiệu quả tính toán và tăng khả năng tổng quát hóa trên dữ liệu mới.

• Thử nghiệm với việc bỏ 4 biến “diastolic” , “systolic” , “gender” và “height_cm”

rf_model1 <- randomForest(class ~ age  + weight_kg + body_fat + grip_force + sit_and_bend_forward_cm + sit_ups_counts+ broad_jump_cm ,
                          data = train_data, ntree = 500, mtry = 3, importance = TRUE)
print(rf_model)

## 
## Call:
##  randomForest(formula = class ~ ., data = train_data, ntree = 500,      mtry = 3, importance = TRUE) 
##                Type of random forest: classification
##                      Number of trees: 500
## No. of variables tried at each split: 3
## 
##         OOB estimate of  error rate: 26.27%
## Confusion matrix:
##      A    B    C    D class.error
## A 1970  282   26    5   0.1371003
## B  525 1371  315   77   0.4007867
## C  199  457 1539  121   0.3354922
## D   27  120  259 1894   0.1765217

• Dự đoán trên tập test

# Dự đoán
rf_predictions <- predict(rf_model, newdata = test_data)

# Đánh giá mô hình
conf_matrix <- table(Predicted = rf_predictions, Actual = test_data$class)
print(conf_matrix)

##          Actual
## Predicted   A   B   C   D
##         A 860 203  81  11
##         B 149 632 176  50
##         C   8 120 643 108
##         D   3  42  67 785

# Tính độ chính xác
accuracy <- sum(diag(conf_matrix)) / sum(conf_matrix)
print(paste("Accuracy:", accuracy))

## [1] "Accuracy: 0.741493143727781"

• Hàm đánh giá tổng hợp

# Gọi hàm đánh giá
results <- eval_multi_class(conf_matrix)

# Hiển thị kết quả đánh giá
print("Evaluation Metrics:")

## [1] "Evaluation Metrics:"

print(results)

## $Precision
##         A         B         C         D 
## 0.8431373 0.6339017 0.6649431 0.8228512 
## 
## $Recall
##         A         B         C         D 
## 0.7445887 0.6276068 0.7315131 0.8751394 
## 
## $F1_Scores
##         A         B         C         D 
## 0.7908046 0.6307385 0.6966414 0.8481902 
## 
## $Macro_Precision
## [1] 0.7412083
## 
## $Macro_Recall
## [1] 0.744712
## 
## $Macro_F1
## [1] 0.7415937
## 
## $Accuracy
## [1] 0.7414931
## 
## $Kappa
## [1] 0.6550024

• Vẽ biểu đồ tương quan cho ma trận nhầm lẫn

library(ggplot2)
library(reshape2)

# Heatmap ma trận nhầm lẫn
conf_matrix <- melt(as.matrix(conf_matrix_qda))
ggplot(conf_matrix, aes(x = Predicted, y = Actual, fill = value)) +
  geom_tile() +
  geom_text(aes(label = value), color = "white") +
  scale_fill_gradient(low = "blue", high = "red") +
  theme_minimal() +
  labs(title = "Confusion Matrix Heatmap", x = "Predicted", y = "Actual")

Dựa vào các hệ số đánh giá mô hình dự đoán thì kết quả sau khi bỏ 4 biến đầu vào không ảnh hưởng đến mô hình. Vì vậy việc bỏ đi các biến trên là hoàn toàn phù hợp.

Kết Luận

Qua việc đánh giá hiệu suất của mô hình và lựa chọn ra các biến có sự ảnh hưởng mạnh mẽ tới biến phân loại class thì ta có thể kết luận rằng:

• 2 biến diastolic và systolic thật sự không ảnh hưởng quá nhiều đến hiệu suất của mô hình, 3 biến height_cm, body_fat và gender có ảnh hưởng một phần nhỏ nhưng không đáng kể nên có thể loại bỏ. Vậy 5 biến này là các biến không ảnh hưởng tới việc phân loại các đối tượng trong class.

• Các biến còn lại như biến sit_ups_counts, broad_jump_cm và age là những biến có ảnh hưởng mạnh mẽ nhất tới việc phân loại.